Mitonuclear and phenotypic discordance in an Atlantic Forest frog hybrid zone

Abstract Discordance between mitochondrial and nuclear DNA is common among animals and can be the result of a number of evolutionary processes, including incomplete lineage sorting and introgression. Particularly relevant in contact zones, mitonuclear discordance is expected because the mitochondrial genome is haploid and primarily uniparentally inherited, whereas nuclear loci are evolving at slower rates. In addition, when closely related taxa come together in hybrid zones, the distribution of diagnostic phenotypic characters and their concordance with the mitochondrial or nuclear lineages can also inform on historical and ongoing dynamics within hybrid zones. Overall, genetic and phenotypic discordances provide evidence for evolutionary divergence and processes that maintain boundaries among sister species or lineages. In this study, we characterized patterns of genetic and phenotypic variation in a contact zone between Cycloramphus dubius and Cycloramphus boraceiensis, two sister species of frogs endemic to the Atlantic Coastal Forest of Brazil. We examined genomic‐scale nuclear diversification across 19 populations, encompassing the two parental forms and a contact zone between them. We compared the distribution of genomic DNA variability with that of a mitochondrial locus (16S) and two morphological traits (dorsal tubercles and body size). Our results reveal multiple divergent lineages with ongoing admixture. We detected discordance in patterns of introgression across the three data types. Cycloramphus dubius males are significantly larger than C. boraceiensis males, and we posit that competition among males in the hybrid zone, coupled with mate choice by females, may be one mechanism leading to patterns of introgression observed between the species.


| INTRODUC TI ON
Animal species often show intraspecific geographic variation associated with diversification (Coyne & Orr, 2004;Endler, 1977;Lee et al., 2019).Even in recently diverged species, microhabitat or topographic variation can isolate lineages, facilitating diversification and speciation (Devitt et al., 2011;Sánchez-Montes et al., 2018).The divergence of populations can lead to premating barriers that prevent interbreeding when these populations are again in sympatry (Pfennig & Rice, 2014;Seddon, 2005).However, these premating isolating barriers are often imperfect and divergent species that subsequently come into secondary contact often interbreed, resulting in a range of possible outcomes (e.g., Maag et al., 2023;Souza et al., 2024).In cases where hybrid offspring occur in zones of secondary contact, postmating barriers may evolve, reinforcing the distinction between species through unfit hybrids, a consequence of parental genomic incompatibility or reduced ecological competitive fitness (Coyne & Orr, 2004;Singhal & Moritz, 2012;Yoshida et al., 2019).Conversely, hybridization and the subsequent backcrossing can promote introgressive hybridization, a process known across a wide range of animal taxa, providing a source of novel and potentially adaptive genetic variation in hybrid populations (Anderson, 1948;Dufresnes & Martínez-Solano, 2020;Grant & Grant, 1994;Lewontin & Birch, 1966;Todesco et al., 2016).Clearly, the outcome of interactions between species across secondary contact zones is highly variable and complex, and understanding these variable outcomes informs on mechanisms leading to the maintenance of divergent lineages or incipient species.
For genetic introgression, distinct selective constraints, such as coding or noncoding loci, effective population sizes, and whether introgression occurs in organellar or autosomal DNA, can lead to genealogical discordances (Maddison, 1997;Singhal & Moritz, 2012).
Discordance between organellar and nuclear loci has been identified in many taxa, particularly in contact zones (Edwards, 2009;Toews & Brelsford, 2012).These discordances are expected because the mitochondrial genome is haploid and typically uniparentally inherited, whereas nuclear loci and phenotypic traits are presumably evolving at slower rates.In some cases, regions of discordance in mitochondrial and nuclear genetic patterns are also coincident with variability in phenotypes that mediate mate selection (Harris et al., 2018;Toews & Brelsford, 2012;Wham et al., 2021).Thus, studies of contact zones have the potential to reveal selective, behavioral, and demographic mechanisms acting on populations, such as differential selection (Cheviron & Brumfield, 2009), sex-biased dispersal (Crochet et al., 2003), mechanisms underlying mate choice (Wirtz, 1999), historical spatial expansion or contraction, and the genetic drift (Currat et al., 2008).
Nonetheless, hybridization between anuran sister species does occur (Haddad et al., 1994).The high diversity in the Atlantic Coastal Forest biome stems from historically dynamic shifts in habitat distributions that resulted in multiple contractions and expansions of species ranges, promoting opportunities for lineage divergences and hybridization among differentiated forms (Amaral et al., 2018;Carnaval et al., 2014;Thomé et al., 2014).Here, we characterized patterns of genetic and phenotypic variation along a contact zone between Cycloramphus dubius and C. boraceiensis, two sister species of frogs endemic to the Atlantic Coastal Forest in southeastern Brazil that have largely parapatric distributions (de Sá, Haddad et al., 2020;Haddad et al., 2013).Our focal species occur in the states of São Paulo and Rio de Janeiro, and have highly similar ecological niches; both are saxicolous, living and breeding in crevices formed by rocks or roots at the margins of fast-flowing streams in forested environments, commonly next to waterfalls (Giaretta & Cardoso, 1995;Giaretta & Facure, 2003;Haddad et al., 2013;Haddad & Prado, 2005;Pedrozo et al., 2024).
Among frogs, females typically choose mates based on territorial and/or phenotypic traits, such as breeding site and/or male body size, therefore enhancing competition among males (de Sá, Consolmagno et al., 2020;Muralidhar et al., 2014;O'Brien et al., 2020;Richardson et al., 2010;Roithmair, 1994).In our two focal species of Cycloramphus, males actively compete and defend territories, where they mate and guard their egg clutches after reproduction, setting the stage for female choice based on these behaviors (de Sá, Haddad et al., 2020;Giaretta & Cardoso, 1995;Giaretta & Facure, 2003;Haddad & Prado, 2005).Cycloramphus dubius and C. boraceiensis are morphologically similar, but diagnosable because C. dubius lacks distinct dorsal white tubercles that are present in C. boraceiensis (Heyer, 1983).Furthermore, both species have similar karyotypes (2n = 26; Beçak et al., 1970;Noleto et al., 2011;Silva et al., 2001) and have partially sympatric distributions, thus providing a good case for the study of gene flow and potential introgression between recently diverged species.
We explore divergence between our two focal frog species with reduced representation genomic data (nuDNA) and compare genomic variation to the distribution of mitochondrial variation (mtDNA), dorsal morphology, and body size along a secondary contact zone.Our goals in this study are to: (1) characterize the contact zone and the outcome of gene flow between our two focal species; and (2) explore potential evolutionary mechanisms that might have governed introgression in different character types.By assessing how secondary contact affects the distribution of genetic and phenotypic diversity, our study advances the understanding of diversification and speciation processes in the Atlantic Forest, a biodiversity hotspot (Marques & Grelle, 2021;Morellato & Haddad, 2000;Myers et al., 2000;Vasconcelos et al., 2014).

| Sample collection
Cycloramphus dubius and C. boraceiensis are saxicolous frogs that are largely parapatric and their ranges come into contact along the central coast of the São Paulo state, near the municipalities of Biritiba Mirim and Salesópolis (de Sá, Haddad et al., 2020).From 2011 to 2016, we sampled populations of both species across their geographic distributions, collecting frogs and tadpoles in the field by hand at night.We supplemented those with tissue samples from previous collections from our laboratories, for a total of 302 individuals from 19 populations.We euthanized field collected frogs and tadpoles with anesthetic overdose (5% lidocaine) and preserved tissues from tadpoles (fins or the whole specimens), juveniles, and adults (livers, toe tips, or thigh muscles) in 100% ethanol.All vouchers and tissues are deposited in the Célio F. B. Haddad amphibian collection (CFBH), Departamento de Biodiversidade, Instituto de Biociências, Universidade Estadual Paulista (UNESP), Rio Claro, São Paulo, Brazil.For these same populations, we collected data on species-level diagnostic dorsal morphologies (for juveniles and adult females and males; described in phenotypic data section below) and body size (for adult females and males).We supplemented phenotypic data from our own field samples with data from museum collections, resulting in data for 721 individuals (Table 1 and Figure 1).Dorsal morphology was characterized for 653 specimens and body size measurements were taken from 371 adults (Table 1).Focusing on the potential contact zone, we sampled outward from there to include most of the areas in which C. dubius and C. boraceiensis occur, only lacking frogs from the extreme northeastern range of C. boraceiensis (in the municipality of Angra dos Reis; Bittencourt -Silva & Silva, 2013).We included samples from the type localities of both species: Alto de Paranapiacaba (site 9, ALTO;
We used the program TCS v. 1.21 (Clement et al., 2000) to generate a haplotype network.The connection limit was set at 95% and gaps in sequences were read as missing data.The haplotype network was visualized using tcsBU (dos Santos et al., 2016).

| Genomic DNA sequencing
Following methods in Peterson et al. (2012), we generated double digest restriction-site associated DNA sequence libraries (ddRADseq) and collected genome-wide single nucleotide polymorphism (SNP) data to characterize nuclear genetic diversification along our hybrid transect (nuDNA).We included uniquely barcoded/ indexed individuals in six libraries (three libraries with 40 individuals each and another three with 56 individuals each), including all those samples for which we had appropriate concentration of DNA.We combined digestion and adapter ligation steps, using 100 ng genomic DNA (10.5 μL at 10 ng μL−1), 15 units of both SbfI-HF and MspI (NEB; 0.75 μL each), 300 units of T4 DNA ligase (NEB, 0.75 μL), 3 μL of 10 × CutSmart ligase buffer (NEB), 1.0 μL each of a uniquely barcoded SbfI-P1 adapter (one of eight) and MspI-P2 adapter (Integrated DNA Technologies; 0.25 μM each), 3 μL of 10 mM ATP (NEB), and 9.25 μL of ddH2O for a total reaction volume of 30 μL.We incubated the restriction-ligation mixture at 37°C for 30 min and then at 20°C for 1 h.We pooled 20 μL from each of the eight barcoded samples.We PCR-amplified each pool using 20 ng P1/P2 adapted DNA template (10 μL at 2 ngμL−1), 12.5 μL 2X Phusion PCR Master Mix (2X, NEB), and 1.25 μL each of forward index primer and a uniquely indexed reverse primer (one of seven; 5 μM each, Integrated DNA Technologies).We pooled PCRs across index groups and size-selected them to retain fragments from 300 to 1,000 bps.Each library was single-end sequenced to 100 bps on one lane of the Illumina HI-SEQ 2500 at the Cornell University Genomics Facility, yielding good quality sequences for 273 individuals (Table 1).

| Quality filtering and SNP calling
Reads were processed in STACKS v. 1.48 (Catchen et al., 2013) with the module process_radtags to discard low-quality reads and reads with ambiguous barcodes or RAD cut sites.The remaining reads were demultiplexed to produce fastq files for each frog individual.
Next, we used the module denovo_map to merge stacks (i.e., identical set of reads) into loci within individuals and to build a catalog of loci across individuals.We set the minimum depth of coverage  (Heyer, 1983).Cycloramphus dubius lacks distinct dorsal white tubercles that are present in C. boraceiensis (white arrow points to a white tubercle).Numbered circles represent sampled populations and are color-coded according to the proportions of the three dorsal morphologies: Orange represents the morphology typical of C. dubius, blue represents the morphology typical of C. boraceiensis, and black represents the intermediate morphology.Numbered circles are sized according to the sample sizes of specimens examined for dorsal morphology (Table 1).Populations 9, 10, and 11 exhibit the three dorsal morphologies, so their circles are enlarged inside the white box for easier visualization.
required to create a stack (−m option) and the maximum number of pairwise differences allowed between any two stacks within a locus (−M option) to three reads.We allowed a maximum of one mismatch when building the catalog of loci across individuals (−n option).Variants were called using the default SNP model with a genotype likelihood ratio test critical value (α) of 0.05.Using the module populations, we filtered SNPs to retain those present in at least two sampling sites (p = 2) and present in 50% of individuals (r = 0.5).We further filtered our data to remove sequence reads with low depth of coverage (<2X per individual) and/or a minor allele frequency less than 0.5%.

| Population genomic analyses: PCAs, structure, hybrid ancestry, and phylogeny
Using R v. 3.5.2(R Core Team, 2018), we conducted a principal component analysis (PCA) on the centered, unscaled nuDNA genotype matrix to visualize and infer patterns of population structure first using the full dataset.Our initial analysis of the full dataset showed clearly that individuals from Ilhabela (site 14, IBEL; a continental island surrounded by sea water, a natural barrier for frogs) are highly isolated and form their own genetic deme.Therefore, we excluded this insular population from some of our downstream analyses, to enhance the visualization of admixture patterns among populations sampled on a transect across the ranges and the hybrid zone of C. dubius and C. boraceiensis.We further quantified population structure and examined nuDNA genetic admixture using Entropy (Gompert et al., 2014), a Bayesian method that deals with uncertainty in nuDNA genotypes and estimates the proportion of ancestry of a set of individuals from contributing populations.Next, we used the adegenet R package v. 2.1.1 to identify the best-fitting number of clusters given the Bayesian Information Criterion (BIC) values for K ranging from 1 to 12.This optimal clustering solution was then used to perform a discriminant analysis of principal components (DAPCs), a multivariate analysis that minimizes within-group variance while maximizing among-group variance (Jombart et al., 2010).Next, we estimated hybrid ancestry between C. dubius and two identified lineages of C. boraceiensis.To determine the ancestry of the putative hybrids, we assigned two source populations and calculated the fraction of nuDNA loci that combine ancestry from these two parental populations in Entropy (Gompert & Buerkle, 2016).We generated triangle plots (Gompert & Buerkle, 2016), where hybrids with different admixture proportions (q) are aligned along the x-axis and ranked according to the fraction of loci at which each individual has ancestry from both parental taxa (Q 12 , interpopulation ancestry, y-axis).This analysis allows us to separate hybrids that are progeny from a cross involving one or both parental taxa (backcrosses or F 1 ) and thus have maximal interpopulation ancestry for a given admixture proportion (on the edges of the triangle; Gompert & Buerkle, 2016), from progeny of crosses between hybrid individuals (F 2 and beyond) that have less than maximal interpopulation ancestry for a given admixture proportion (middle of the triangle; Gompert & Buerkle, 2016).The insular population (site 14, IBEL) was excluded from the hybrid ancestry analysis.
To understand the phylogenetic relationships between C. dubius and C. boraceiensis populations, we first built a maximum likelihood (ML) phylogeny using RAxML-NG v. 1.2.2 (Kozlov et al., 2019) in the CIPRES Science Gateway v. 3.3 (Miller et al., 2010).We used the full SNP dataset, including the insular population.The best-fit nucleotide substitution model for our dataset was determined by the corrected Akaike Information Criterion (AICc) in the jModelTest v. 2.1.8.(Darriba et al., 2012).We applied the GTR + I + G nucleotide substitution model, with 20 starting trees (10 random and 10 parsimony trees), and assessed branch support by running 1000 bootstrap replicates.We used FigTree v. 1.4.4 to edit the final tree (Rambaut, 2018), apllying midpoint rooting.
We also inferred population trees using a coalescent model and Bayesian MCMC approach implemented in SNAPP (Bryant et al., 2012).SNAPP is a tool developed for phylogenetic analyses based on biallelic markers (such as SNPs), directly inferring species trees by integrating over all possible genealogies rather than sampling them.This approach provides high statistical power but is computationally expensive, with complexity scaling according to the number of samples and markers (Yoder et al., 2013).For our analyses, we used SNAPP as implemented in BEAST v. ), for a total of 48 individuals (Rodger et al., 2021;Stetter & Schmid, 2017).We used the preliminary ML topology (Figure S1) to guide our sample choice and also included in the reduced dataset those samples that were not grouping with their own populations.This necessarily biases topologies and parameter estimates in our cloudogram but also represents the full discordance among genealogies of hybridizing lineages.
We assessed the output using Tracer v. 1.7.2 (Rambaut et al., 2018) to ensure acceptable mixing and convergence, confirming effective sample sizes >200 for all parameters.The two runs for each replicate were combined using LogCombiner v. 2.7.4 (Drummond & Rambaut, 2007), using a burn-in of 10%.The posterior distribution of gene (SNP) trees was visualized using DensiTree v. 2.6.1 (also implemented in BEAST; Bouckaert, 2010) to capture uncertainty in topology and branch lengths.Finally, we summarized the data in a maximum clade credibility (MCC) tree using TreeAnnotator v. 2.7.4 (also implemented in BEAST; Drummond & Rambaut, 2007), after discarding the first 20% of samples as burn-in.
de SÁ et al.

| Phenotypic data and analyses
At each sampling site, we collected data on dorsal morphology and adult body size (snout-vent length, SVL).We augmented our phenotypic dataset by scoring dorsal morphologies and measur-  1).We also discovered a novel combination of the two basic dorsal morphology types, which we categorized as "intermediate", in which the dorsum is partially covered by white tubercles.We calculated the proportion of individuals of the three dorsal types at each sampling site.
We measured adult SVL with calipers (±0.1 mm).To estimate sexual size dimorphism, we used a size dimorphism index, calculated as the ratio of female to male body size subtracted by one [SDI: (female SVL/male SVL) - 1; Lovich & Gibbons, 1992].Our data were normally distributed (Shapiro-Wilk test), and showed equality of variances (Levene's test), so we compared averages of female and male body sizes, and sexual size dimorphism among genetic lineages using t-test (t) or one-way ANOVA test (F) followed by Tukey test (Q'), in R v. 3.5.2(R Core Team, 2018).We considered statistical significance when p < .05(Zar, 1998).

| Geographic cline analyses
To

| Genomic structure and admixture
Our final nuDNA genetic data set included 5972 SNPs with an average of 153× depth of coverage per individual per SNP (±51 SD) and 14.5% missing data.The genetic clustering found in the PCA (Figure 3a) mirrored the geographic distribution of populations (Figure 1).Specifically, PC1, which explained 20.2% of the as the optimal number of clusters found within the data (Figure S3).
Additional visual inspection of K values ranging from 7 to 8 revealed no additional population clusters, instead highlighting the extent of admixture among populations (Figure S2).

| Hybrid ancestry
We generated three pairwise triangle plots of hybrid ancestry for the three lineages represented among continental populations (Figure 4).We excluded the isolated insular deme on Ilhabela because it showed little evidence of admixture with neighboring demes.Thus, our triangle plots represent hybrids between C. du-  Buerkle, 2016).Numbered populations are defined in Table 1.

| Phylogenetic relationships
The  5a).This result is confirmed by the MCC analysis (Figure 5b).

| Dorsal morphologies and body sizes
Populations varied in the proportions of the three dorsal morphology types (Figure 1).Populations are fully defined in Table 1.
de SÁ et al.
morphologies.We found that females are undifferentiated by their sizes (t = 1.84, p = .06869),but that males with the typical C. dubius dorsal morphology are significantly larger than those with a C. boraceiensis morphology (t = 2.16, p = .03526;Table 2 and Figure 6).
Commonly, when mainland populations have geographic structure Note: Values are shown as average (in mm) ±SD; in parentheses range and sample size (n).We report Q'-values and p-values from statistical comparisons (Tukey test) of female-female and male-male sizes between different lineages.We also show size variations focusing on the hybrid zone (sites 1-13) but grouping adult females and males by their dorsal morphologies (i.e., independently from their DNA assignments) and also excluding intermediate morphologies.Numbered populations are defined in Table 1.
and are large, and/or when island isolation happened recently, the small population isolated on the island renders the mainland populations paraphyletic (Avise, 2000;Funk & Omland, 2003;Johnson et al., 2000;Moritz, 1994).Most continental islands off the Atlantic Coastal Forest were connected to the mainland for the last time during the last glacial maximum, ~13 thousand years ago (Fleming et al., 1998;Souza et al., 2005;Suguio & Martin, 1978).Over multiple periods of sea-level oscillations and marine introgressions, these islands became fully isolated during the Holocene, when the relative sea level rose from approximately −60 m to levels that we observe today (Fleming et al., 1998;Souza et al., 2005;Suguio & Martin, 1978).These sea level fluctuations promoted different levels of isolation across taxa, with some organisms reaching speciation via vicariance (Brasileiro et al., 2007;Grazziotin et al., 2006;Marques et al., 2002).
Geotectonic structures might explain phylogeographic breaks commonly observed in the Atlantic Forest (Amaral et al., 2013;Amaro et al., 2012) as Cycloramphus, also show a split in the region of São Sebastião (Duryea et al., 2015;Sabbag et al., 2018).Likewise, the terrestrially breeding Ischnocnema parva frog complex from the family Brachycephalidae (Gehara et al., 2017) and the aquatically breeding Paratelmatobius poecilogaster frog complex from the family Leptodactylidae (Santos et al., 2020) have lineages splitting around São Sebastião.As we accumulate evidence from phylogeographic datasets for different taxa, we will illuminate general mechanisms leading to the diversification in this biodiversity hotspot.

| Taxonomic remarks
Our two focal species have high levels of similarity in morphology, behaviors, and ecological niches, and are also closely related (de Sá, Haddad et al., 2020;Giaretta & Cardoso, 1995;Giaretta & Facure, 2003;Heyer, 1983;Pedrozo et al., 2024).Nonetheless, they are phenotypically distinguishable (Heyer, 1983) and our mitochondrial data confirm that they belong to independent haplogroups, even with partial sympatry.Heyer (1983) wrote that "C.(Frost, 2023;Heyer, 1983), and the holotypes show the set of diagnostic phenotypes for each species.Our results uncover an interesting taxonomic challenge, because both type localities fall within the hybrid zone we characterized in this study.We did not have access to DNA samples from the type locality of C. dubius and the two individuals of C. boraceiensis we had from its type locality were genomically categorized as hybrids.Today, C. boraceiensis is extinct at its type locality (Lopes et al., 2020), thus, we do not know if specimens in the type series are hybrids or not.Overall, the fact that our hybrid zone is relatively narrow and that parental forms diagnosed by morphology are found even in populations within the hybrid zone led us to conclude that these two species are maintaining distinct identities, despite genetic introgression.A future re-evaluation focused on the taxonomy of these sister species is clearly necessary.

| Insights on evolutionary mechanisms operating in this discordant hybrid zone
We examined variation in dorsal morphology for C. dubius and C. boraceiensis, and detected more diversity than has been previously described (Heyer, 1983) (Funk & Omland, 2003).
Hybridization is not uncommon among closely related frog species and has been observed even in species that show elaborate courtship behaviors or variation in sexually selected phenotypic traits, conditions that likely mediate assortative mating (Akopyan et al., 2020;Nali et al., 2022).In the Cycloramphus hybrid zone described here, populations of the two species do not appear to face geographic, ecological, behavioral, or biological barriers to reproduction.Our findings suggest that pre-or post-zygotic reproductive isolation mechanisms are failing to prevent hybridization between these two sister species in secondary contact.
We did detect subtle but statistically significant differences in body sizes between the species.Notably, females are the same size, but males are variable, which may indicate distinct selective forces acting on each sex (de Sá, Haddad et al., 2020).First, larger males found on Ilhabela (site 14, IBEL) might be a consequence of insular effects (Stamps & Buechner, 1985).Second, in the hybrid zone, males of C. dubius and males of southwestern C. boraceiensis (C.boraceiensis L1) are similar in sizes; however, when examining only males with species-typical dorsal morphologies and grouping these males only by those morphologies, we found that males with typical C. dubius dorsal morphologies are significantly larger than those with C. boraceiensis dorsal morphologies.We know that both C. dubius and C. boraceiensis have similar reproductive strategies and behaviors, with males aggressively defending territories that encompass the best egg-laying sites and guarding their eggs (Giaretta & Cardoso, 1995;Giaretta & Facure, 2003), and with females likely choosing mates based on male traits and territories (de Sá, Haddad et al., 2020;Lipshutz, 2018;O'Brien et al., 2020;Santana et al., 2020;Valencia-Aguilar et al., 2020).The main hypothesis we posit here is that sexual selection might be an evolutionary factor contributing to the distribution of hybrids in this system.Based on the body size variation of males in the hybrid zone, it is possible that, if any interspecific male-male competition takes place, C. dubius would have an advantage, benefiting from their larger size when defending resources and excluding competing males (see de Sá, Consolmagno et al., 2020;de Sá, Haddad et al., 2020;Lipshutz, 2018;Muralidhar et al., 2014;Pearson & Rohwer, 2000), and potentially displacing males of C. boraceiensis.Two admixing species differing in body sizes could lead to non-reciprocal hybridization.Males of the larger species might be chosen by females of the smaller species, but males of the smaller species might be rejected by females of the larger species (Grant & Grant, 1997;Lipshutz, 2018;Wirtz, 1999).Larger males of C. dubius might be acting as supernormal stimulus for both conspecific and C. boraceiensis females, whereas smaller males of C. boraceiensis might be acting as a subnormal stimulus for females (de Sá, Haddad et al., 2020;Lipshutz, 2018;Richardson et al., 2010;Wells, 2010;While et al., 2015;Winkelmann et al., 2014).If females do prefer larger males, female choice would also favor nonreciprocal introgression.
Alternatively, factors other than the male body size or advertisement calls can potentially influence the introgression dynamics detected in this hybrid zone.For example, genetic compatibilities, with females preferring allospecific males because of heterosis benefits or other genetic-related advantages (e.g., Marr et al., 2002) or chemical communication, with females preferring allospecific males due to chemical signals (e.g., Johansson & Jones, 2007).
Intrasexual competition and differential male fitness is a wellknown selective mechanism in the genus Cycloramphus and in other species in the family (de Sá, Consolmagno et al., 2020;de Sá, Haddad et al., 2020;Muralidhar et al., 2014).Our previous work on the genus Cycloramphus demonstrated that competition among males is pervasive and an evolutionary driver of male body size (de Sá, Haddad et al. 2020).Here, we propose that male-male competition coupled with female mate choice may also play an important role in shaping the distribution of organellar and genomic variation, and phenotypic variation in a secondary contact zone between two divergent sister species.Future studies examining hybrid zones should also consider the role that sexual selection mechanisms might play in shaping anuran diversification.
2.7.6 (Bouckaert et al., 2019) via CIPRES.We used the default model and prior parameters (models: μ = 1, v = 1, coalescence rate = 10; priors: α = 11.75, β = 109.73,κ = 1, and λ = 0.01) and ran two independent replicate MCMC runs of 1000,000 iterations with sampling every 1000 steps and a burn-in of 10%.To ensure feasible computational times, we reduced the SNP dataset to 2408 SNPs, randomly selected from the full dataset (of 5972 SNPs), and included only 12 randomly chosen individuals from each lineage (C.dubius, C. boraceiensis L1, insular C. boraceiensis, and C. boraceiensis L2 estimate the extent of introgression between C. dubius and C. boraceiensis, we estimated clinal transitions of mitochondrial haplotype frequencies, allele frequencies and admixture proportions from the SNP dataset, and dorsal morphology frequencies.We fit genetic and phenotypic data to geographic cline models using the Metropolis-Hasting Markov chain Monte Carlo algorithm implemented in the R package hzar(Derryberry et al., 2014).We excluded the insular population (site 14, IBEL) from the geographic cline analyses.To generate clines, we calculated pairwise Euclidean distances (in km) from the southern-most sampling site from geographic coordinates with the program Geographic Distance Matrix Generator v. 1.2.3(Ersts, 2011).Mitochondrial haplotype data were coded as 1 or 0, with 1 indicating that the haplotype belonged to the C. dubius mitochondrial clade and 0 indicating that the haplotype belonged to one of the two mitochondrial lineages of C. boraceiensis.Clinal transitions of admixture proportions from the SNP dataset (q) were estimated based on the C. dubius deme, so individuals with high q had high C. dubius ancestry, and individuals with low q had high C. boraceiensis ancestry.Therefore, q values indicate proportions of an individual's genome attributed to each source population(Shastry et al., 2021).Additionally, we estimated clinal transitions of dorsal morphologybased on the proportion of individuals with the uniform dorsal morphology typical of C. dubius.We fit clines with a null model and 15 cline models based on all possible combinations of trait interval (fixed to 0 and 1, observed values or estimated values) and tail fits (none fitted, left only, right only, mirrored, or both estimated separately).Model selection was performed using corrected Akaike Information Criterion (AICc) scores and the ML model parameters (cline center and width) were extracted from the best model.Coincidence of cline centers and concordance of cline widths were determined using confidence intervals of two log-likelihood scores.

|
Mitochondrial haplotype networkWe identified a total of 43 mtDNA haplotypes, including 11 haplotypes of the C. dubius lineage and 32 of two C. boraceiensis lineages (Figure 2).Haplotypes of C. dubius and C. boraceiensis differed by a minimum of nine mutational steps.The 11 haplotypes of C. dubius are present at sites 1-12.The 32 haplotypes of C. boraceiensis consisted of 19 haplotypes present at sites 12-14 (hereafter called lineage Cycloramphus boraceiensis L1, which includes the insular population) that differed by a minimum of five mutational steps from the 13 haplotypes at sites 15-19 (hereafter called lineage Cycloramphus boraceiensis L2).The haplotype network shows one of the C. boraceiensis L1 haplotypes shared with C. dubius at one sampling locality (site 12, SALE).The lineage C. boraceiensis L1 is more genetically distinct from C. dubius haplotypes than the lineage C. boraceiensis L2 (Figure 2); however, C. boraceiensis L1 is geographically closer to the contact zone between the two species (Figure 1).
Figure 3c), corresponding to the four clusters revealed by the PCA analysis.Further partitioning the data into six populations (K = 6, Figure S2) revealed additional structure by distinguishing the sampling sites at the extreme ends of our sampling transect.Based on the K values evaluated using BIC scores, DAPC identified K = 6-10

Figure 4 ).
Figure 4).We found little evidence of late generation hybrids, which fall within the center of the triangle, having intermediate levels of genome-wide ancestry but lower interpopulation ancestry (Figure 4).Most hybrids are backcrossed individuals between C. dubius and C. boraceiensis L1, and between C. dubius and C. boraceiensis L2, as indicated by the points falling toward the left and right sides of the triangles (Figure4).We found evidence of more de SÁ et al. backcrossed hybrids with a pure C. dubius parent than any other deme.Finally, individuals that fell along the x-axis (intermediate admixture proportions, with no interpopulation ancestry) are interpreted as showing a pattern of geographic isolation by distance (see Gompert & Buerkle, 2016).
inferred topology from the SNP-based RAxML-NG analysis recovers major clades in agreement with our other analyses, with separation between C. dubius and C. boraceiensis, and population structure within C. boraceiensis (C.boraceiensis L1, insular C. boraceiensis, and C. boraceiensis L2).However, the overall topology has relatively low support, especially at the base of the tree (Figure S1).The insular C. boraceiensis population is clearly isolated with strong nodal support (site 14, IBEL).Some individuals from site 1 (PTOL; sample 226), site 8 (CUBA; samples 138 and 158), and site 12 (SALE; samples 271 and 273) fall at the base of the C. dubius clade and may be the result of admixture or incomplete lineage sorting between C. dubius and C. boraceiensis L1.Similarly, individuals from site 1 (PTOL; samples 222 and 225) and site 3 (ITAN; sample 268) are potentially admixed between C. dubius and C. boraceiensis L2, whereas individuals from site 13 (SSEB; sample 201) and site 15 (CARA; sample 248) are potentially admixed between C. boraceiensis L1 and L2.Our cloudogram generated via DensiTree analysis, using a reduced SNP dataset and with only 48 individuals, shows the overall disagreement between trees in the SNAPP run.The only two lineages with strong support are the insular C. boraceiensis and C. boraceiensis L2.All other lineages show topological uncertainty in their relationships.Cycloramphus dubius, C. boraceiensis L1, and C. boraceiensis L2 form a three-way polytomy at the base of the tree, with some topologies supporting closer relationships between C. dubius and C. boraceiensis L1 (Figure Figure 6).Lastly, focusing particularly on the hybrid zone detected between C. dubius and C. boraceiensis L1 (sites 1-13), we also grouped adult females and adult males only based on their dorsal morphologies (i.e., independently from their DNA assignments), also excluding intermediate F I G U R E 5 Populations trees for Cycloramphus dubius and Cycloramphus boraceiensis (C.boraceiensis L1, C. boraceiensis insular population, and C. boraceiensis L2) using a reduced SNP dataset and with only 48 individuals.(a) The posterior distribution of gene trees indicate overall ambiguous relationships among populations, as showed in our cloudogram.(b) SNAPP provides full posterior support only for the insular C. boraceiensis and C. boraceiensis L2, as showed our MCC tree.
phology and mitochondrial haplotype frequencies, and a free trait interval based on observed values was selected for q.The center of the mitochondrial haplotype cline was estimated near SALE (site 12; Figure1) at 158.1 km (CI: 147.8-161.8)from the western-most site, with a width of 0.8 km (CI: 0.2-17.1).For q, the cline center was estimated near BIRI (site 11; Figure1) at 133.8 km from the western-most site, with a width of 28.2 km (CI: 0.6-65.6).For dorsal morphology, the cline center was estimated farther west at 107.7 km (CI: 103.1-112.4),near ALTO (site 9; Figure1), with a width of 46.6 km (CI: 38.2-57.2).The log-likelihood values for the cline centers do not overlap.The center of the mitochondrial cline was estimated near the species boundaries, between BIRI (site 11) and SALE (site 12), whereas the q cline center was displaced nearly 24-km west, indicating mitonuclear discordance in the contact zone between Cycloramphus boraceiensis and C. dubius.The dorsal morphology cline was displaced yet nearly 26-km west of the q cline center, indicating that the white-tubercled dorsal morphology typical of C. boraceiensis is also discordant from the previous two cline centers (Figure 7).evidence of phenotypic and genetic variation among populations of C. dubius and C. boraceiensis, two sister species of frogs endemic to the Atlantic Coastal Forest in southeastern Brazil.MtDNA and nuDNA patterns for C. dubius are more homogeneous across populations (sites 1-11) and are recovered as a single cluster.Conversely, C. boraceiensis populations (sites 12-19) are split along their geographic distribution into three main nuDNA clusters: southwestern (lineage C. boraceiensis L1), insular (Ilhabela; site 14, IBEL), and northeastern (lineage C. boraceiensis L2) lineages.The four recovered clusters are further supported by the ddRADseq-based PCAs and the phylogenetic analyses, but with important evidence of introgression among lineages.The insular mtDNA haplotypes from Ilhabela (site 14) fall within the southwestern cluster (C.boraceiensis L1 haplogroup) with some haplotype sharing with coastal populations.In contrast, the nuDNA data show that admixture between insular and C. boraceiensis L1 populations is minor.This pattern of insular genetic differentiation is typical of a recent separation, likely with insular C. boraceiensis populations isolated by sea-level oscillations dur-

boraceiensis
most closely resembles C. dubius", but also indicated that C. dubius lacks distinct dorsal, white-tipped tubercles.We verified this taxonomic diagnostic characters by examining the dorsal morphology of seven C. dubius and 29 C. boraceiensis individuals from the type series (housed at the MZUSP) as well as samples from outside the contact zone, and concluded that the two species are in fact diagnosable throughout much of their ranges based on dorsal tubercles, but in addition they also have differences in other phenotypic characters, with C. dubius showing a more granular dorsal texture and with C. boraceiensis showing slightly more developed foot webbing.The designated type locality for C. dubius is Alto da Serra, Paranapiacaba, Santo André, São Paulo, and the type locality for C. boraceiensis is Estação Biológica de Boracéia, Salesópolis, São Paulo Few advertisement calls have been characterized for C. dubius and C. boraceiensis across their distributions, but the few calls that have been recorded show some interesting differences.Calls of the two species have a similar range of dominant frequencies, but C. dubius males from Santos (site 7; SANT) emit longer calls (0.2 s) than C. boraceiensis L1 males from Salesópolis (site 12; SALE) or C. boraceiensis L2 from Ubatuba (site 16; UBAT), which emit advertisement calls of 0.03 and 0.06 s in length, respectively
Note: We report the locality (population) and final sample sizes for genetic and phenotypic data included in downstream analyses.Populations are numbered according to the map in Figure1.Abbreviations: RJ, Rio de Janeiro state; SP, São Paulo state.F I

G U R E 1
Population sampling of Cycloramphus dubius and Cycloramphus boraceiensis in the states of São Paulo and Rio de Janeiro, southeastern Brazil.Based on mtDNA and nuDNA, C. dubius occurs at sites 1-11, lineage C. boraceiensis L1 at sites 12-13, insular C. boraceiensis population at site 14, and C. boraceiensis L2 at sites 15-19.Numbered populations are defined in Table 1.Images are representatives of each species, which are distinguished based on dorsal morphology sizes of specimens at four institutions in Brazil: CFBH, Rio Size variations of adult females and males from each genetic lineage, Cycloramphus dubius (sites 1-11), southwestern lineage Cycloramphus boraceiensis L1 (sites 12-13), insular C. boraceiensis population (site 14), and northeastern lineage C. boraceiensis L2 (sites 15-19).
TA B L E 2 . Overall, the distributions of C. dubius and C. Body sizes by sex and SDIs (=size dimorphism indices; see 2. Materials and Methods for details), all grouped by recovered genetic clusters.Cycloramphus dubius occurs from PTOL to BIRI (sites 1-11); southwestern Cycloramphus boraceiensis lineage (C.boraceiensis L1) occurs in SALE and SSEB (sites 12-13); insular C. boraceiensis population occurs in IBEL (site 14); and northeastern C. boraceiensis lineage (C.boraceiensis L2) occurs from CARA to PARA (sites 15-19).See Figure 1 for site location on map.Numbered populations and sample sizes are fully defined in Table 1.Whereas females are all the same size, males show more variability and statistical differences (different sizes separated by letters A, B, and C on boxes; see results for details).SÁ et al. both species have their geographic limits that closely match hydrographic basins of the Brazilian coast, with C. dubius located in the Santos Lowlands basin and C. boraceiensis in the North Coast basin.The topography of the Serra do Mar might also explain diversification within C. boraceiensis, because Serra do Mar runs parallel to the Brazilian coastline, including high escarpments closer to the ocean near the city of São Sebastião, São Paulo state (site 13, SSEB; boraceiensis along the Brazilian Atlantic Coastal Forest are associated with escarpments of the Serra do Mar Mountain chain, and F I G U R E 6 F I G U R E 7 Best-supported cline models for mitochondrial haplotypes (in red; mtDNA haplotype), genome-wide admixture proportions based on K = 1 (in yellow; SNP admixture, q), and (in blue) dorsal morphologies of Cycloramphus dubius and Cycloramphus boraceiensis.de . Except for a few individuals with intermediate morphologies, most individuals in the three C. boraceiensis genetic northeastern-southwestern nuDNA cline between C. dubius and C. boraceiensis, with introgression and hybrid individuals (sites 9-12).The nuDNA-based PCAs and hybrid ancestry plots further confirm introgression between the two species.In addition, Salesópolis (site 12, SALE) is the only sampled population with C. dubius and C. boraceiensis (C.boraceiensis L1) individuals sharing mitochondrial